Brownian motion and Diffusion

import numpy as np
import scipy 

from numpy.random import normal, choice, uniform

import ipywidgets as widgets
import matplotlib.pyplot as plt
plt.style.use('fast')

%matplotlib inline

Mean square displacement (MSD) of a random walker

• After time n number of steps (or time t) how far has random walker moved from the origin?

\[Z_n = \sum^{n-1}_{i=0}X_n\]

• We quantify this by computing Mean Square Displacement (MSD). Note that the mean is computed over N number of simulated trajectories (ensemble average). Invoking central limit theorem, or simply realizing that off diagonal terms drop off we end up with expected 1/2 saling with steps.

\[ MSD(n)= \Big\langle \big ( Z_n - Z_0 \big)^2 \Big \rangle \sim n^{1/2} \]

Mean square displacement (MSD) and Diffusion

• Einstein developed a theory of diffusion based on random walk ideas and obtained a key equation relating mean square displacement of particle in d dimensions to time \((t = n \Delta t)\).

• Consider recording n number of positions for a particle in solution starting from origin \(r_0=0\)

\[R_n = \sum^{n}_{i=0} r_i\]

• Repeating experiment N number of times we can compute ensemble average

\[\langle R^2_n\rangle = \sum_i \sum_j \langle r_i r_j \rangle = \sum_i \langle r^2_i \rangle = \sum_i d \cdot \langle \delta x^2_i \rangle = d \cdot n \cdot \delta x^2\]

• Where we decomposed displacement in d independent dimensions.

• We next write number of steps in terms of time increments \(n=t/\delta t\)

\[\langle R^2_n\rangle = d \cdot \frac{t}{\delta t} \cdot \delta x^2\]

• Grouping constants together we defined the diffusion coefficient:

\[\langle R^2 (t) \rangle = 2d \cdot D \cdot t\]

\[D = \frac{\langle \delta x^2 \rangle}{2\delta t}\]

• We end up with a general expression for a mean square displacement as a function of time. Any motion which adheres to this scaling will be called “diffusive”

\[MSD(t) = 2d D \cdot t^{1/2}\]

MSD of a 1D random walk

def rw_1d(T, N):
    '''
    L: trajectory length
    N: Number of trajecotries
    returns np.array with shape (T, N) 
    '''
    
    # Create random walks 
    r  = choice([-1,1], size=(T, N))
    
    #Accumulate position
    rw = r.cumsum(axis=0)

    #Set initial position 
    rw[0,:]=0 
    
    return rw

T, N = 2000, 1000
rw = rw_1d(T, N)

t = np.arange(T)

R2 = (rw[:, :] - rw[0, :])**2 # Notice we subtract initial time

msd =  np.mean(R2, axis=1)    # Notice we average over N

plt.loglog(t, np.sqrt(msd), lw=3) 

plt.loglog(t, np.sqrt(t), '--')

plt.title('Compute mean square deviation of 1D random walker',fontsize=15)
plt.xlabel('Number of steps, n',fontsize=15)
plt.ylabel(r'$MSD(n)$',fontsize=15);

MSD of a 2D random walk

def rw_2d(T, N):
    '''2d random walk function:
    T: trajectory length
    N: Number of trajecotry
    returns np.array with shape (T, N)
    '''
    verteces = np.array([(1,  0),
                         (0,  1),
                         (-1, 0),
                         (0, -1)])
    
    rw       = verteces[choice([0,1,2,3], size=(T, N))]
    
    rw[0, :, :] = 0
    
    return rw.cumsum(axis=0)

traj = rw_2d(T=10000, N=1000)

#Simulate 2D random walk
T, N = 10000, 1000
traj = rw_2d(T, N)

#Compute RSD 
dx = (traj[:, :, 0]- traj[0, :, 0]) 
dy = (traj[:, :, 1]- traj[0, :, 1]) 

R2     = np.mean(dx**2 + dy**2, axis = 1)   # notice how we averaging over N

fig, ax  = plt.subplots(nrows=2, figsize=(10,10))

t = np.arange(T) # time axis
ax[1].loglog(t, np.sqrt(R2), lw=3, alpha=0.5);
ax[1].loglog(t, np.sqrt(t), '--');

ax[0].set_title('2D random walker',fontsize=15)
ax[0].plot(traj[:3000, :5, 0], traj[:3000, :5, 1]);
ax[0].set_xlabel('X')
ax[1].set_xlabel('Y')

ax[1].set_xlabel('Number of steps, n',fontsize=15)
ax[1].set_ylabel(r'$MSD(n)$',fontsize=15);

Brownian motion

• The Brownian motion describes the movement of a particle suspended in a fluid resulting from random collisions with the quick molecules in the fluid (diffusion).

• A small particle undergoes large number of molecular collisions when going from one step to another or after \(dt\) time. As a result each displacement over time \(dt\) can be viewed as a sum of ranomd collisions which can be approximated by a normal distribution via Central Limit Theorem.

• Thus to simulate brownian motion we draw ranom displacements from normally distribution.

\[x(t+dt)-x(t)=N(0,\sqrt{2D dt})\]

We assume we have started at position \(\mu=0\) and our variance is given by \(\sigma^2=2Dt\) Where D is the diffusion coefficnets which is related to parameters of discree random walk as shown in the lecture.

\[x(t+dt)=x(t)+\sqrt{2D dt} \cdot N(0,1)\]

In the last step we re-wrote brownian motion equation in a convenient way by shifting normally distributed radnom variable by \(\mu\) and scaling it by \(\sigma\)

\[N(\mu, \sigma^2) = \mu + \sigma N(0,1) \]

def brown(T, N, dt=1, D=1):
    
    """
    Creates 3D brownian path given:
    time T 
    N=1 trajecotires
    dt=1 timestep
    D=1 diffusion coeff
    returns np.array with shape (N, T, 3)
    """
    
    n = int(T/dt) # how many points to sample
    
    dR = np.sqrt(2*D*dt) * np.random.randn(N, n, 3) # 3D position of brownian particle
    
    R = np.cumsum(dR, axis=1) # accumulated 3D position of brownian particle
    
    return R

• Below we proceed to simulate continuous random walk in 1D-3D

• We will plot trajecotires and distributions of brownian particle using interactive plotting via ipywidgets and holoviews/plotly interface.

R=brown(T=3000, N=1000)
print(R.shape)

(1000, 3000, 3)

@widgets.interact(t=(10,3000-1))
def brownian_plot(t=10):
    
    fig, ax = plt.subplots(ncols=2)
    
    ax[0].plot(R[:5, :t, 0].T, R[:5, :t, 1].T);
    
    ax[1].hist(R[:, 10, 0], density=True, color='red')
    ax[1].hist(R[:, t, 0], density=True)
    
    ax[1].set_ylim([0,0.1])
    
    ax[0].set_ylim([-200, 200])
    ax[0].set_xlim([-200, 200])
    
    fig.tight_layout()

import holoviews as hv
hv.extension('plotly')

plots = []
for i in range(10):
    
    plot = hv.Path3D(R[i,:,:], label='3D random walk').opts(width=600, height=600, line_width=5)
    plots.append(plot)
    
hv.Overlay(plots) 

rw_curve = [hv.Curve((R[i,:,0], R[i,:,1])) for i in range(10)]

xdist = hv.Distribution(R[:,10,0], ['X'], ['P(X)'])
ydist = hv.Distribution(R[:,10,1], ['Y'], ['P(Y)'])

hv.Overlay(rw_curve) << ydist << xdist

Diffusion Equation

The movement of individual random walkers \(\leftrightarrow\) density of walkers \(\rho(\vec{r},t)\)

Diffusion equation:

Formulated empirically as Fick’s laws

\[\frac{\partial\rho}{\partial t} = \mathcal{D}\nabla^2\rho\]

• This is a 2nd order PDE! Unlike equations of motion diff eq shows irreersibile behaviour

• This one exactly solvable. But in general reaction-diffusion PDEs difficult to solve analytically.

• Can solve numerically by writing derivatives as finite differences!

• Can also simulate via random walk!

• Diffusion coefficient \(D\), Units \([L^2]/[T]\)

Important special case solution (here written in 1d):

\[\rho(x,t) = \frac{1}{\sqrt{2\pi \sigma(t)^2}}\exp\left(-\frac{x^2}{2\sigma(t)^2}\right),\]

where \(\sigma(t)=\sqrt{2{D}t}\)

• density remains Gaussian

• Gaussian becomes wider with time

• check that this is indeed a solution by plugging into the diffusion equation!

def sigma(t, D = 1):
    return np.sqrt(2*D*t)

def gaussian(x, t):
    return  1/np.sqrt(2*np.pi*sigma(t)**2) * np.exp(-x**2/(2*sigma(t)**2)) #

@widgets.interact(t=(1,100,1))
def diffusion(t=0.001):
    
    R = brown(T=101, N=1000)
    x = np.linspace(-20, 20, 100)
    
    plt.plot(x, gaussian(x, 1), '--', color='orange', label='t=0')
    
    plt.plot(x, gaussian(x, t), lw=3, color='green', label=f't={t}')
    
    plt.hist(R[:,t,0], density=True, alpha=0.6, label='simulation hist')
    
    plt.legend()
    plt.ylabel('$p(x)$')
    plt.xlabel('$x$')
    plt.xlim([-25, 25])

References

The mighty little books

• “Random Walks in Biology”, H Berg (1993)

• “Physical models of Living systems”, P Nelson (2015)

More in depth

• “Simple Brownian Diffusion: An Introduction to the Standard Theoretical Models”, D Gillespie

• “Stochastic Processes for Physicists” K Jacobs

On the applied side

• Brownian Motion: Elements of Colloid Dynamics A P Philipse (2018)